// --*- C++ -*------x-----------------------------------------------------------
//
//
// Class:           LogAverage
//
// Author:          Enrico Negri
//
// Project name:    BiocompUP Align
//
// Date:            12/2007
//
// Description:     Calculate scores for profile to profile alignment using
//                  sum of pairs method. Some explanations can be found in:
//
//                  Guoli Wang, Roland L. Dunbrack jr.
//                  Scoring profile-to-profile sequence alignments.
//                  Institute for Cancer Research, Fox Chase Cancer Center,
//                  Philadelphia, Pennsylvania 19111, USA. March 16, 2004.
//
// -----------------x-----------------------------------------------------------

              /*************************************************
               * sub is BLOSUM62 matrix standard log-odds form *
               *************************************************/

#include <LogAverage.h>

namespace Biopool
{

// CONSTRUCTORS:

LogAverage::LogAverage(SubMatrix *sub, Profile *pro1, Profile *pro2)
	: ScoringFunction(), sub(sub), pro1(pro1), pro2(pro2)
{ }


LogAverage::LogAverage(const LogAverage &orig) : ScoringFunction(orig)
{
	copy(orig);
}


LogAverage::~LogAverage()
{ }


// OPERATORS:

LogAverage&
LogAverage::operator = (const LogAverage &orig)
{
	if (&orig != this)
		copy(orig);
	return *this;
}


// PREDICATES:

double
LogAverage::scoringSeq(int i, int j)
{
	double s = 0.00;

	for (AminoAcidCode amino1 = ALA; amino1 <= TYR; amino1++)
	{
		double tmp = 0.00;

		for (AminoAcidCode amino2 = ALA; amino2 <= TYR; amino2++)
			tmp += exp(sub->score[aminoAcidOneLetterTranslator(amino1)]
				[aminoAcidOneLetterTranslator(amino2)]) *
					pro2->getAminoFrequencyFromCode(amino2, (j-1));

		s += (pro1->getAminoFrequencyFromCode(amino1, (i-1)) * tmp);
	}

	return log(s);
}


// MODIFIERS:

void
LogAverage::copy(const LogAverage &orig)
{
	ScoringFunction::copy(orig);
	sub = orig.sub->newCopy();
	pro1 = orig.pro1->newCopy();
	pro2 = orig.pro2->newCopy();
}


LogAverage*
LogAverage::newCopy()
{
	LogAverage *tmp = new LogAverage(*this);
	return tmp;
}

} // namespace
